Notebook arguments

• sigma (float): standard deviation of additive Gaussian noise to be simulated
• time_window (float): seconds, integration window duration
• time_step (float): seconds, time step for the moving integration window
• time_start (float): seconds, start of time axis (kinetics starts at t = t0).
• time_stop (float): seconds, stop of time axis (kinetics starts at t = t0).
• t0_vary (bool): whether models should vary the curve origin (t0) during the fit
• true_params (dict): parameters used to generate simulated kinetic curves
• num_sim_cycles (int): number of times fit is repeated (Monte-Carlo)
• taus (tuple): list of values for time-costant tau simulated during repeated fits (Monte-Carlo).

---

Simulated Kinetic Curve Fit

This notebook fits simulated exponential transients with additive Gaissian noise in order to study time-constant fitting accuracy. In particular we compare a simple exponential model with a more realistic model with integration window, checking the effect on the fit results.

You can either run this notebook directly, or run it through the master notebook for batch processing.

Imports

Fitting models

Models used to fit the data.

1. Simple Exponential

In this model, we define the model function as an exponential transient:

$$ y = f(t) = A \cdot e^{-t/\tau} + K$$

The python function implementing it is:

• models.exp_func().

Next cell defines and initializes the fitting model (lmfit.model.Model) including parameters' constrains:

2. Integrated Exponential

A more realistic model needs to take into account that each data point is the result of an integration over a time window $w$:

$$f(t) = A \cdot e^{-t/\tau} + K$$$$y(t) = \int_{t}^{t+w} f(t')\;dt'$$

In other words, when we process a measurement in time chunks, we are integrating a non-stationary signal $f(t)$ over a time window $w$. This integration causes a smoothing of $f(t)$, regardless of the fact that time is binned or is swiped-through with a moving windows (overlapping chunks).

Numerically, $t$ is discretized with step equal to (time_step / decimation).

The python function implementing this model function is:

• models.expwindec_func().

And, finally, we define and initialize the fitting model parameters' constrains:

Generative model

These are the models used to generate the simulates (noisy) data.

1. Simple Exponential + Noise

In this simple model, we simulate random data $Y$ as an exponential decay plus additive Gaussian noise:

$$ Y(t_k) = f(t_k) + N_k $$$$ \{N_k\} \sim {\rm Normal}\{\mu=0; \sigma\}$$$$ \Delta t = t_k - t_{k-1} = \texttt{time_step}$$

2. Integrated Exponential + Noise

For the "integrating window" model, we first define a finer time axis $\theta_i$ which oversamples $t_k$ by a factor $n$. Then we define the function $Y_f$ adding Gaussian noise $\sqrt{n}\,N_i$, with $n$ times larger variance:

$$ Y_f(\theta_i) = f(\theta_i) + \sqrt{n}\,N_i $$$$ \Delta \theta = \theta_i - \theta_{i-1} = \texttt{time_step} \;/\; n$$

Finally, by averaging each time window, we compute the data on the coarse time axis $t_k$:

$$ Y_w(t_k) = \frac{1}{m}\sum_{i} Y_f(\theta_i)$$

Here, for each $t_k$, we compute the mean of $m$ consecutive $Y_f$ values. The number $m$ is chosen so that $m\, \Delta \theta$ is equal to the time window.

Noise amplitude

The amplitude of the additive noise ($\sigma$) is estimated from the experimental kinetic curves. In particular we take the variance from the POST period (i.e. the steady state period after the transient). The POST period has been chosen because it exhibits higher variance than the PRE period (i.e. the steady state period before the transient). These values have been calculated in 8-spot bubble-bubble kinetics - Summary.

In both models we define the noise amplitude as sigma (see first cell):

sigma = 0.016

Time axis

We also define the parameters for the time axis $t$:

time_start = -900    # seconds 
time_stop  =  900    # seconds
time_step  = 5       # seconds

Kinetic curve paramenters

The simulated kinetic curve has the following parameters:

true_params = dict(
    tau = 60,          # time constant
    init_value = 0.3,  # initial value (for t < t0) 
    final_value = 0.8, # final value   (for t -> +inf)
    t0 = 0)            # time origin

Single kinetic curve fit

Here we simulate one kinetic curve and fit it with the two models (simple exponential and integrated exponential).

Draw simulated data

Time axis for simulated data:

An ideal transient (no noise, no integration):

A simulated transient (including noise + integration):

Plot the computed curves:

Fit data

Fit the "Integrated Exponential" model:

Fit the "Simple Exponential" model:

Print and plot fit results:

Monte-Carlo Simulation

Here, fixed the model paramenters, we generate and fit several noisy datasets. Then, by plotting the distribution of the fitted parameters, we assess the stability and accuracy of the fit.

Parameters

The number simulation cycles is defined by num_sim_cycles. Current value is:

The fixed kinetic curve parameters are:

While tau is varied, taking the following values:

Functions

Here we define two functions:

• draw_samples_and_fit() draws a set of data and fits it with both models

• monte_carlo_sim() run the Monte-Carlo simulation: calls draw_samples_and_fit() many times.

NOTE: Global variables are used by previous functions.

Run Monte-Carlo simulation

Run the Monte-Carlo fit for a set of different time-constants (taus) and save results in two DataFrames, one for each model.

Results1 - Simple Exponential

Results2 - Integrated Exponential

Conclusions

The last two multipanel figures compare the fitting accuracy of the model parameter for the simple-exponential and integrated-exponential models. We note that, in particular for the tau parameter, the integrated exponential model is significantly more accurate, providing good estimates at much smaller integration times.

This comparison demonstrates empirically the strong advantage in using the theoretically more correct integrated exponential model.